표준화된 두 변수 간의 공변하는 관계를 나타내는 통계량이다. 일반적으로는 영국의 통계학자 피어슨(K. Pearson, 1857-1936)이 개발한 피어슨 적률상관계수(product moment correlation)를 가리키는 명칭이다. 이것은 선형적 관계를 나타내며 -1과 +1 사이에 있고 아래의 공식으로 계산된다.
상관계수는 상관되는 변수들의 성질에 따라 종류가 다양하다. Pearson의 적률상관계수는 두변수가 수량변수(연속변수)일 때 사용한다. 그러나 두 변수 모두 자연스런 이분변수일 때의 파이계수(phi coefficient), 두 변수가 모두 자연스런 다분변수일 때의 스피어만 서열상관(Spearman rank order correlation), 하나의 수량변수와 다른 하나의 자연스런 이분변수 간의 점이연상관(point-biserial correlation), 수량변수와 자연스런 다분변수 간의 점다연상관(point-polyserial correlation)은 모두 피어슨 적률상관계수의 공식을 간편화시켜서 계산하므로 피어슨상관계열(Pearsonian correlations)이라 한다.
피어슨 상관계열이 아닌 특수상관은 -1과 +1의 범위를 벗어나며, 사분상관(tetrachoric correlation), 이연상관(biserial correlation), 다연상관(poly-serial correlation) 등이 있다.
상관계수의 종류 참고 사이트 link
install.packages(‘ggplot2’) install.packages(‘Hmisc’)
install.packages(‘ggm’)
install.packages(‘polycor’)
install.packages(‘boot’)
img <- list()
img <- paste0("image/f16-",c(1:14),".jpg")
imglist <- NULL
for (i in c(1:14)) imglist <- paste0(imglist, " ", "", sep="")
library(ggplot2) #ggplot()
library(ggm) #rcorr() 함수 있으나 지저분한 결과 나온다고 함. pcor()
library(Hmisc) #rcorr() 이 library 추천하고 있음.
library(polycor)
library(boot)
RBGL이 없다고 뜨는 경우 ###source(“http://bioconductor.org/biocLite.R”); biocLite(c(“graph”,“RBGL”,“Rgraphviz”))
pearson: 연속 또는 등간및 비율척도 상관계수
spearman or kendall: 순위상관계수
cor()는 상관계수만 보여주고, rcorr()는 상관계수와 함께 n과 p-value 보여주나 kendall안됨.
####cor.test()는 모든 걸 보여주고 모든게 됨. p216 Table 6.2 참조.
E = read.delim("6correlation data/Exam Anxiety.dat", header = T)
head(E)
## Code Revise Exam Anxiety Gender
## 1 1 4 40 86.298 Male
## 2 2 11 65 88.716 Female
## 3 3 27 80 70.178 Male
## 4 4 53 80 61.312 Male
## 5 5 4 40 89.522 Male
## 6 6 22 70 60.506 Female
###cor(E, use = 'complete.obs', method = 'pearson') #error
cor(E$Exam, E$Anxiety, use = 'complete.obs', method = 'pearson')
## [1] -0.4409934
cor(E$Exam, E$Anxiety, use = 'complete.obs', method = 'kendall')
## [1] -0.2847919
cor(E$Exam, E$Anxiety, use = 'pairwise.complete.obs', method = 'kendall')
## [1] -0.2847919
rcorr(E$Exam, E$Anxiety, type = 'pearson')
## x y
## x 1.00 -0.44
## y -0.44 1.00
##
## n= 103
##
##
## P
## x y
## x 0
## y 0
###rcorr(E, type = 'pearson') #error
cor.test(E$Exam, E$Anxiety, method = 'pearson')
##
## Pearson's product-moment correlation
##
## data: E$Exam and E$Anxiety
## t = -4.938, df = 101, p-value = 3.128e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.5846244 -0.2705591
## sample estimates:
## cor
## -0.4409934
cor.test(E$Exam, E$Anxiety, alternative = 'less', method = 'pearson')
##
## Pearson's product-moment correlation
##
## data: E$Exam and E$Anxiety
## t = -4.938, df = 101, p-value = 1.564e-06
## alternative hypothesis: true correlation is less than 0
## 95 percent confidence interval:
## -1.0000000 -0.2995071
## sample estimates:
## cor
## -0.4409934
cor.test(E$Exam, E$Anxiety, alternative = 'less', method = 'pearson', conf.level = .99)
##
## Pearson's product-moment correlation
##
## data: E$Exam and E$Anxiety
## t = -4.938, df = 101, p-value = 1.564e-06
## alternative hypothesis: true correlation is less than 0
## 99 percent confidence interval:
## -1.0000000 -0.2362782
## sample estimates:
## cor
## -0.4409934
이 파트는 세 변수 이상의 상관계수를 공부.
rcorr()는 matrix형태(수치변수만) 받음.
E2 = E[,c('Revise', 'Exam', 'Anxiety')]
cor(E2)
## Revise Exam Anxiety
## Revise 1.0000000 0.3967207 -0.7092493
## Exam 0.3967207 1.0000000 -0.4409934
## Anxiety -0.7092493 -0.4409934 1.0000000
EMatrix = as.matrix(E[,c('Revise', 'Exam', 'Anxiety')])
rcorr(EMatrix)
## Revise Exam Anxiety
## Revise 1.00 0.40 -0.71
## Exam 0.40 1.00 -0.44
## Anxiety -0.71 -0.44 1.00
##
## n= 103
##
##
## P
## Revise Exam Anxiety
## Revise 0 0
## Exam 0 0
## Anxiety 0 0
#cor.test(E$Exam, E$Anxiety) 중복된 것. rcorr(변수 세 개)로 가능.
#cor.test(E$Revise, E$Anxiety) 중복된 것. rcorr(변수 세 개)로 가능.
cor(E2)^2 #상관계수 제곱한 것.
## Revise Exam Anxiety
## Revise 1.0000000 0.1573873 0.5030345
## Exam 0.1573873 1.0000000 0.1944752
## Anxiety 0.5030345 0.1944752 1.0000000
비모수(모집단 분포 신경 안 씀) 상관분석
거짓말 대회 등수(position)와 그 사람의 창의력(설문조사함, creativity)
liarData = read.delim('6correlation data/The Biggest Liar.dat', header = T)
head(liarData)
## Creativity Position Novice
## 1 53 1 0
## 2 36 3 1
## 3 31 4 0
## 4 43 2 0
## 5 30 4 1
## 6 41 1 0
cor(liarData$Position, liarData$Creativity, method = 'spearman')
## [1] -0.3732184
liarMatrix = as.matrix(liarData[,c('Position','Creativity')])
rcorr(liarMatrix)
## Position Creativity
## Position 1.00 -0.31
## Creativity -0.31 1.00
##
## n= 68
##
##
## P
## Position Creativity
## Position 0.0111
## Creativity 0.0111
#단측검정: 대립가설(H1, alternative)을 rho<0 이다로 설정하고 검정. 귀무가설은 rho>=0 인데, p<0.05이므로 귀무가설을 기각함.
cor.test(liarData$Position, liarData$Creativity, method = 'spearman', alternative ='less')
## Warning in cor.test.default(liarData$Position, liarData$Creativity, method
## = "spearman", : Cannot compute exact p-value with ties
##
## Spearman's rank correlation rho
##
## data: liarData$Position and liarData$Creativity
## S = 71948, p-value = 0.0008602
## alternative hypothesis: true rho is less than 0
## sample estimates:
## rho
## -0.3732184
#창의력(creativity)높아서 거짓말을 잘할수록 position(1이 1등임)은 낮아져야 하니까
cor(liarData$Position, liarData$Creativity, method = 'kendall')
## [1] -0.3002413
cor.test(liarData$Position, liarData$Creativity, method = 'kendall', alternative = 'less')
##
## Kendall's rank correlation tau
##
## data: liarData$Position and liarData$Creativity
## z = -3.2252, p-value = 0.0006294
## alternative hypothesis: true tau is less than 0
## sample estimates:
## tau
## -0.3002413
부트스트랩(bootstrap) 또는 부트스트래핑은 “현재 상황에서 어떻게든 한다”는 뜻이다. 또, 사물의 초기 단계에서 단순 요소로부터 복잡한 체계를 구축하는 과정을 가리키는 경우도 있다.
[네이버 지식백과] 부트스트랩법 [bootstrap method, ~法] (생명과학대사전, 초판 2008., 개정판 2014., 도서출판 여초)
자료를 통해 얻은 통계량의 표본오차를 확률분포의 가정을 두지 않고 논파라메트릭적으로 평가하기 위한 방법의 하나. 잭나이프법이나 무작위화 검정과 같이 컴퓨터 사용을 전제로 한 새로운 통계방법의 하나이다. 주어진 데이터 세트를 원래의 모집단을 대표하는 독립표본으로 가정하고, 그 자료로부터의 중복을 허용한 무작위 재추출을 컴퓨터로 복수의 자료를 작성하여 각각에서 얻은 통계량(부트스트랩반복)을 계산한다.
그리고 얻어진 복수의 통계량의 오차분산을 구한다. 부트스트랩법은 생물학에서의 사용도 차차 확산되어 가고 있지만 특히 계통추정론 분야에서는 계통수의 신뢰성을 평가할 목적으로 널리 사용하고 있다. 개개의 형질을 자료점이라고 생각하고, 관찰된 형질집합을 1,000회 내지 10,000회 정도 부트스트랩하여 각각의 형질집합으로부터 계통수를 추정한다. 얻어진 복수의 계통수 전체에서 가지(단계통군)의 출현율(부트스트랩확률)을 계산한다. 출현율이 높은 가지(가령 90% 이상)일수록 신뢰할 수 있는 단계통군을 반영한다고 볼 수 있다.
#bootstrapping correlation, spearman과 pearson에 대해서도 적용 가능.
#i는 i번째 bootstrapping set을 의미함.
bootTau <- function(liarData, i){
cor(liarData$Position[i], liarData$Creativity[i], method = 'kendall', use = 'complete.obs')
}
boot_kendall <- boot(liarData, bootTau, 2000) #object 생성: boot(data, function, replications)
boot_kendall
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = liarData, statistic = bootTau, R = 2000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* -0.3002413 2.045563e-05 0.09350031
boot.ci(boot_kendall) #95% confidence interval
## Warning in boot.ci(boot_kendall): bootstrap variances needed for
## studentized intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_kendall)
##
## Intervals :
## Level Normal Basic
## 95% (-0.4835, -0.1170 ) (-0.4890, -0.1193 )
##
## Level Percentile BCa
## 95% (-0.4812, -0.1115 ) (-0.4727, -0.0955 )
## Calculations and Intervals on Original Scale
boot.ci(boot_kendall, conf = 0.99)
## Warning in boot.ci(boot_kendall, conf = 0.99): bootstrap variances needed
## for studentized intervals
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_kendall, conf = 0.99)
##
## Intervals :
## Level Normal Basic
## 99% (-0.5411, -0.0594 ) (-0.5479, -0.0604 )
##
## Level Percentile BCa
## 99% (-0.5401, -0.0526 ) (-0.5165, -0.0199 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstable
한 쪽 변수가 양적인 변수(구간, 비율)이고, 다른 쪽 변수는 이진형(binary, dichotomous)일 경우, 두 변수의 관계를 측정할 때 사용된다.
catData = read.csv('6correlation data/pbcorr.csv', header = T)
head(catData)
## time gender recode
## 1 41 1 0
## 2 40 0 1
## 3 40 1 0
## 4 38 1 0
## 5 34 1 0
## 6 46 0 1
cor.test(catData$time, catData$gender)
##
## Pearson's product-moment correlation
##
## data: catData$time and catData$gender
## t = 3.1138, df = 58, p-value = 0.002868
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.137769 0.576936
## sample estimates:
## cor
## 0.3784542
cor.test(catData$time, catData$recode)
##
## Pearson's product-moment correlation
##
## data: catData$time and catData$recode
## t = -3.1138, df = 58, p-value = 0.002868
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.576936 -0.137769
## sample estimates:
## cor
## -0.3784542
###polyserial
catFrequencies = table(catData$gender)
prop.table(catFrequencies)
##
## 0 1
## 0.5333333 0.4666667
polyserial(catData$time, catData$recode)
## [1] -0.4749256
#구문: pcor(c("var1", "var2", "control1", "control2" etc.), var(dataframe))
pc = pcor(c('Exam', 'Anxiety', 'Revise'), var(E2))
pc
## [1] -0.2466658
pc^2
## [1] 0.06084403
#구문: pcor.test(pcor object, number of control variables, sample size)
pcor.test(pc, 1, 103)
## $tval
## [1] -2.545307
##
## $df
## [1] 100
##
## $pvalue
## [1] 0.01244581
link 준부분상관계수는 예측하고자 하는 준거변인(Y)을 설명하는 예언변인(X1,X2)의 고유한 설명량을 확인할 수 있는 것이다. 이 책의 hierarchical regression, analysis of covariance을 공부하면 이 내용을 다시 확인할 수 있다. 여기서는 이 정도만 하겠다.
교과서에 없고, 홈페이지의 Oliver Twisted를 참고하면 나온다. 아래가 그 사이트의 pdf문서임. page 7과 page 10 부근을 참조하면 된다.
maleExam<-subset(E, Gender == "Male", select= c("Exam", "Anxiety"))
femaleExam<-subset(E, Gender == "Female", select= c("Exam", "Anxiety"))
cor(maleExam)
## Exam Anxiety
## Exam 1.0000000 -0.5056874
## Anxiety -0.5056874 1.0000000
cor(femaleExam)
## Exam Anxiety
## Exam 1.0000000 -0.3813845
## Anxiety -0.3813845 1.0000000
zdifference<-function(r1, r2, n1, n2){
zd<-(atanh(r1)-atanh(r2))/sqrt(1/(n1-3)+1/(n2-3)) #atanh: arc tangent h(쌍곡삼각함수)
p <-1 - pnorm(abs(zd))
print(paste("Z Difference: ", zd))
print(paste("One-Tailed P-Value: ", p))
}
zdifference(-0.506, -0.381, 52, 51)
## [1] "Z Difference: -0.768709306290097"
## [1] "One-Tailed P-Value: 0.221032949510287"
r’=0.5*ln(abs((1+r)/(1-r))), z=(r1’-r2’)/sqrt((1/(N1-3))+1/(N2-3)) link
tdifference<-function(rxy, rxz, rzy, n){
df<-n-3
td<-(rxy-rzy)*sqrt((df*(1 + rxz))/(2*(1-rxy^2-rxz^2-rzy^2+(2*rxy*rxz*rzy))))
p <-pt(td, df)
print(paste("t Difference: ", td))
print(paste("One-Tailed P-Value: ", p))
}
tdifference(-0.441, -0.709, 0.397, 103)
## [1] "t Difference: -5.09576822523987"
## [1] "One-Tailed P-Value: 8.21913727738007e-07"
r값(correlation coefficient) 자체가 Effect size를 의미한다.
주의: Caution for non-parametric statistic
Spearman은 순위(rank)를 공유하는 변화량의 정도를 보여주는 것. 정규분포(Normal distribution)에 가까울수록 \(r^2\)은 pearson만큼 의미가 있다.
Kendall tau는 두 변수 간 분산의 공유 정도를 나타내지 않는다. 제곱하여 \(tau^2\) 형식으로 사용할 수 없다. 통상 평가자들 간의 점수의 일치도를 의미한다.
Kendall tau는 66~75%정도 작기 때문에 Pearson이나 Spearman과 비교할 수 없다.
point-biserial/biserial 판단에 따라서 상관정도가 달라질 수 있으므로 변수의 특성을 명확히 파악하여야 한다.